#ifdef MPIF

MODULE MPIComm_Mod
!USE MPI
IMPLICIT NONE
include "mpif.h"
!USE UtilFunction,     ONLY : ArrayLoc
INTEGER :: MPI_MC_COMM, MPI_CMFD_COMM
INTEGER, PARAMETER, PRIVATE :: nmax = 5000000
INTEGER :: ierr

INTERFACE BCAST
  MODULE PROCEDURE BCAST_REAL0D
  MODULE PROCEDURE BCAST_REAL1D
  MODULE PROCEDURE BCAST_REAL2D
  MODULE PROCEDURE BCAST_REAL3D
  MODULE PROCEDURE BCAST_LOG0D
  MODULE PROCEDURE BCAST_INT0D
  MODULE PROCEDURE BCAST_INT1D
  MODULE PROCEDURE BCAST_INT2D
  MODULE PROCEDURE BCAST_INT3D
END INTERFACE

INTERFACE REDUCE
  MODULE PROCEDURE REDUCE_REAL0D
  MODULE PROCEDURE REDUCE_REAL1D
  MODULE PROCEDURE REDUCE_REAL2D
  MODULE PROCEDURE REDUCE_REAL3D
  MODULE PROCEDURE REDUCE_REAL4D
  MODULE PROCEDURE REDUCE_REAL5D
  MODULE PROCEDURE REDUCE_INT0D
  MODULE PROCEDURE REDUCE_INT1D
  MODULE PROCEDURE REDUCE_INT2D
  MODULE PROCEDURE REDUCE_INT3D
END INTERFACE

INTERFACE REDUCEnBCAST
  MODULE PROCEDURE REDUCEnBCAST_Real0D
  MODULE PROCEDURE REDUCEnBCAST_Int0D
END INTERFACE

CONTAINS
FUNCTION nCommStep(n)
INTEGER :: n, nCommStep
INTEGER :: i, j
i = n / nmax
j = mod(n, nmax)
IF(J .EQ. 0) THEN
 nCommStep = i
ELSE
 nCommStep = i + 1
ENDIF
END FUNCTION

SUBROUTINE ArrayLoc(idx, n, Dims, Locs)
IMPLICIT NONE
INTEGER :: n, idx
INTEGER :: Dims(n), LOCS(n)
INTEGER :: m, i, j, k, nn
m = 1
DO i = 1, n
  m = m * Dims(i) 
ENDDO
nn = idx
DO i = n, 1, -1 
  m = m / Dims(i)
  IF(i .EQ. 1) THEN
    LOCS(1) = nn
    EXIT
  ENDIF
  DO j = 1, DIMS(i)
    IF(nn .LE. j * m) EXIT
  ENDDO
  j = nn / m;   k = MOD(nn, m)
  IF(k .EQ. 0) THEN
    j = j
  ELSE
    j = j + 1
  ENDIF  
  nn  = nn - m * (j-1)
  LOCS(i) = j 
  IF(nn .eq. 0 .AND. i .NE. 1) THEN
    LOCS(1:i-1) = DIMS(1:i-1)
    EXIT
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_REAL5D(Dat, DatBuf, n1, n2, n3, n4, n5, comm, lall)
IMPLICIT NONE
REAL(8) :: DAT(n1, n2, n3, n4, n5)
REAL(8) :: DatBuf(n1, n2, n3, n4, n5)
INTEGER :: n1, n2, n3, n4, n5, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ie3, ie4, ie5, ib1, ib2, ib3, ib4, ib5
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(5), Locs(5)
LOGICAL :: lAll
ndat = n1 * n2 * n3 * n4 * n5
m = nCommStep(ndat)
Dims = (/ n1, n2, n3, n4, n5 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 5, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2); ib3 = Locs(3); ib4 = Locs(4); ib5 = Locs(5)
!  CALL ArrayLoc(ie, 3, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2); ie3 = Locs(3)
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(Dat(ib1, ib2, ib3, ib4, ib5), DatBuf(ib1, ib2, ib3, ib4, ib5), n, &
                    MPI_DOUBLE_PRECISION, MPI_SUM, 0, Comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(Dat(ib1, ib2, ib3, ib4, ib5), DatBuf(ib1, ib2, ib3, ib4, ib5), n, &
                      MPI_DOUBLE_PRECISION, MPI_SUM, Comm, ierr)
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_REAL4D(Dat, DatBuf, n1, n2, n3, n4, comm, lall)
IMPLICIT NONE
REAL(8) :: DAT(n1, n2, n3, n4)
REAL(8) :: DatBuf(n1, n2, n3, n4)
INTEGER :: n1, n2, n3, n4, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ie3, ie4, ib1, ib2, ib3, ib4
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(4), Locs(4)
LOGICAL :: lAll
ndat = n1 * n2 * n3 * n4
m = nCommStep(ndat)
Dims = (/ n1, n2, n3, n4 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 4, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2); ib3 = Locs(3); ib4 = Locs(4)
!  CALL ArrayLoc(ie, 3, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2); ie3 = Locs(3)
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(Dat(ib1, ib2, ib3, ib4), DatBuf(ib1, ib2, ib3, ib4), n, &
                    MPI_DOUBLE_PRECISION, MPI_SUM, 0, Comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(Dat(ib1, ib2, ib3, ib4), DatBuf(ib1, ib2, ib3, ib4), n, &
                      MPI_DOUBLE_PRECISION, MPI_SUM, Comm, ierr)
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_REAL3D(Dat, DatBuf, n1, n2, n3, comm, lall)
IMPLICIT NONE
REAL(8) :: DAT(n1, n2, n3)
REAL(8) :: DatBuf(n1, n2, n3)
INTEGER :: n1, n2, n3, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ie3, ib1, ib2, ib3
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(3), Locs(3)
LOGICAL :: lAll
ndat = n1 * n2 * n3
m = nCommStep(ndat)
Dims = (/ n1, n2, n3 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 3, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2); ib3 = Locs(3) 
!  CALL ArrayLoc(ie, 3, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2); ie3 = Locs(3)
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(Dat(ib1, ib2, ib3), DatBuf(ib1, ib2, ib3), n, &
                    MPI_DOUBLE_PRECISION, MPI_SUM, 0, Comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(Dat(ib1, ib2, ib3), DatBuf(ib1, ib2, ib3), n, &
                      MPI_DOUBLE_PRECISION, MPI_SUM, Comm, ierr)
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_INT3D(DAT, DatBuf, n1, n2, n3, comm, lAll)
IMPLICIT NONE
INTEGER :: DAT(n1, n2, n3)
INTEGER :: DatBuf(n1, n2, n3)
INTEGER :: n1, n2, n3, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ie3, ib1, ib2, ib3
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(3), Locs(3)
LOGICAL :: lAll 
ndat = n1 * n2 * n3
m = nCommStep(ndat)
Dims = (/ n1, n2, n3 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 3, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2); ib3 = Locs(3) 
!  CALL ArrayLoc(ie, 3, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2); ie3 = Locs(3)
  IF(.NOT. lALl) THEN
    CALL MPI_REDUCE(Dat(ib1, ib2, ib3), DatBuf(ib1, ib2, ib3), n, &
                    MPI_INTEGER, MPI_SUM, 0, Comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(Dat(ib1, ib2, ib3), DatBuf(ib1, ib2, ib3), n, &
                    MPI_INTEGER, MPI_SUM, Comm, ierr)  
  ENDIF
ENDDO
END SUBROUTINE


SUBROUTINE REDUCE_REAL2D(DAT, DatBuf, n1, n2, comm, lAll)
IMPLICIT NONE
REAL(8) :: Dat(n1, n2)
REAL(8) :: DatBuf(n1, n2)
INTEGER :: n1, n2, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ib1, ib2
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(2), Locs(2)
LOGICAL :: lAll
ndat = n1 * n2
m = nCommStep(ndat)
Dims = (/ n1, n2 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 2, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2)
!  CALL ArrayLoc(ie, 2, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2)
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(Dat(ib1, ib2), DatBuf(ib1, ib2), n, MPI_DOUBLE_PRECISION, &
                    MPI_SUM, 0, Comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(Dat(ib1, ib2), DatBuf(ib1, ib2), n, MPI_DOUBLE_PRECISION, &
                    MPI_SUM, Comm, ierr)  
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_INT2D(DAT, DatBuf, n1, n2, comm, lAll)
IMPLICIT NONE
INTEGER :: DAT(n1, n2)
INTEGER :: DatBuf(n1, n2)
INTEGER :: n1, n2, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ib1, ib2
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(2), Locs(2)
LOGICAL :: lAll
ndat = n1 * n2
m = nCommStep(ndat)
Dims = (/ n1, n2 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 2, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2)
!  CALL ArrayLoc(ie, 2, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2)
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(Dat(ib1, ib2), DatBuf(ib1, ib2), n, &
                    MPI_INTEGER, MPI_SUM, 0, Comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(Dat(ib1, ib2), DatBuf(ib1, ib2), n, &
                    MPI_INTEGER, MPI_SUM, Comm, ierr)  
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_REAL1D(Dat, DatBuf, ndat, comm, lAll)
REAL(8) :: DAT(ndat)
REAL(8) :: DATBuf(ndat)
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n, m
LOGICAL :: lAll
m = nCommStep(ndat)

DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(dat(ib), DatBuf(ib), n, MPI_DOUBLE_PRECISION, MPI_SUM, 0, comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(dat(ib), DatBuf(ib), n, MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
  ENDIF
  
ENDDO

END SUBROUTINE

SUBROUTINE REDUCE_INT1D(Dat, DatBuf, ndat, comm, lAll)
INTEGER :: DAT(ndat)
INTEGER :: DatBuf(ndat)
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n, m
LOGICAL :: lAll
m = nCommStep(ndat)

DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  IF(.NOT. lAll) THEN
    CALL MPI_REDUCE(dat(ib), DatBuf(ib), n, MPI_INTEGER, MPI_SUM, 0, comm, ierr)
  ELSE
    CALL MPI_ALLREDUCE(dat(ib), DatBuf(ib), n, MPI_INTEGER, MPI_SUM, comm, ierr)
  ENDIF
ENDDO
END SUBROUTINE

SUBROUTINE REDUCE_REAL0D(Dat, DatBuf, comm, lAll)
REAL(8) :: DAT
REAL(8) :: DATBuf
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n
LOGICAL :: lAll
n = 1
IF(.NOT. lAll) THEN
  CALL MPI_REDUCE(dat, DatBuf, n, MPI_DOUBLE_PRECISION, MPI_SUM, 0, comm, ierr)
ELSE
  CALL MPI_ALLREDUCE(dat, DatBuf, n, MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
ENDIF


END SUBROUTINE

SUBROUTINE REDUCE_INT0D(Dat, DatBuf, comm, lAll)
INTEGER :: DAT
INTEGER :: DatBuf
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n
LOGICAL :: lAll
n = 1
IF(.NOT. lAll) THEN
  CALL MPI_REDUCE(dat, DatBuf, n, MPI_INTEGER, MPI_SUM, 0, comm, ierr)
ELSE
  CALL MPI_ALLREDUCE(dat, DatBuf, n, MPI_INTEGER, MPI_SUM, comm, ierr)
ENDIF
END SUBROUTINE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE BCAST_REAL3D(DAT, n1, n2, n3, comm)
IMPLICIT NONE
REAL(8) :: DAT(n1, n2, n3)
INTEGER :: n1, n2, n3, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ie3, ib1, ib2, ib3
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(3), Locs(3)
ndat = n1 * n2 * n3
m = nCommStep(ndat)
Dims = (/ n1, n2, n3 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 3, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2); ib3 = Locs(3) 
!  CALL ArrayLoc(ie, 3, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2); ie3 = Locs(3)
  CALL MPI_BCAST(Dat(ib1, ib2, ib3), n, MPI_DOUBLE_PRECISION, 0, Comm, ierr)
ENDDO
END SUBROUTINE

SUBROUTINE BCAST_INT3D(DAT, n1, n2, n3, comm)
IMPLICIT NONE
INTEGER :: DAT(n1, n2, n3)
INTEGER :: n1, n2, n3, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ie3, ib1, ib2, ib3
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(3), Locs(3)
ndat = n1 * n2 * n3
m = nCommStep(ndat)
Dims = (/ n1, n2, n3 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 3, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2); ib3 = Locs(3) 
!  CALL ArrayLoc(ie, 3, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2); ie3 = Locs(3)
  CALL MPI_BCAST(Dat(ib1, ib2, ib3), n, MPI_INTEGER, 0, Comm, ierr)
ENDDO
END SUBROUTINE


SUBROUTINE BCAST_REAL2D(DAT, n1, n2, comm)
IMPLICIT NONE
REAL(8) :: DAT(n1, n2)
INTEGER :: n1, n2, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ib1, ib2
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(2), Locs(2)
ndat = n1 * n2
m = nCommStep(ndat)
Dims = (/ n1, n2 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 2, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2)
!  CALL ArrayLoc(ie, 2, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2)
  CALL MPI_BCAST(Dat(ib1, ib2), n, MPI_DOUBLE_PRECISION, 0, Comm, ierr)
ENDDO
END SUBROUTINE

SUBROUTINE BCAST_INT2D(DAT, n1, n2, comm)
IMPLICIT NONE
INTEGER :: DAT(n1, n2)
INTEGER :: n1, n2, ndat, n, m
INTEGER :: ib, ie, ie1, ie2, ib1, ib2
INTEGER :: comm
INTEGER :: i
INTEGER :: Dims(2), Locs(2)
ndat = n1 * n2
m = nCommStep(ndat)
Dims = (/ n1, n2 /)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL ArrayLoc(ib, 2, Dims, Locs)
  ib1 = Locs(1); ib2 = Locs(2)
!  CALL ArrayLoc(ie, 2, Dims, Locs)
!  ie1 = Locs(1); ie2 = Locs(2)
  CALL MPI_BCAST(Dat(ib1, ib2), n, MPI_INTEGER, 0, Comm, ierr)
ENDDO
END SUBROUTINE

SUBROUTINE BCAST_REAL1D(Dat, ndat, comm)
REAL(8) :: DAT(ndat)
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n, m
m = nCommStep(ndat)
DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL MPI_BCAST(dat(ib), n, MPI_DOUBLE_PRECISION, 0, comm, ierr)
ENDDO

END SUBROUTINE

SUBROUTINE BCAST_INT1D(Dat, ndat, comm)
INTEGER :: DAT(ndat)
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n, m
m = nCommStep(ndat)

DO i = 1, m
  ib = nmax * (i - 1) + 1
  ie = nmax * i; ie = min(ie, ndat)
  n = ie - ib + 1
  CALL MPI_BCAST(dat(ib), n, MPI_INTEGER, 0, comm, ierr)
ENDDO
END SUBROUTINE

SUBROUTINE BCAST_REAL0D(Dat, comm)
REAL(8) :: DAT
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n
ndat = 1
n=1
CALL MPI_BCAST(dat, n, MPI_DOUBLE_PRECISION, 0, comm, ierr)

END SUBROUTINE

SUBROUTINE BCAST_LOG0D(Dat, comm)
LOGICAL :: DAT
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n
ndat = 1
n=1
CALL MPI_BCAST(dat, n, MPI_LOGICAL, 0, comm, ierr)

END SUBROUTINE

SUBROUTINE BCAST_INT0D(Dat, comm)
INTEGER :: DAT
INTEGER :: ndat, comm
INTEGER :: i, j, ib, ie
INTEGER :: n
ndat = 1
n=1
CALL MPI_BCAST(dat, n, MPI_INTEGER, 0, comm, ierr)

END SUBROUTINE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE GetMPIFile(Dat, Master, comm)
IMPLICIT NONE
CHARACTER(80) :: Dat 
LOGICAL :: Master
INTEGER:: comm, ierr
CALL MPI_BCAST(DAT, 80, MPI_CHARACTER, 0, comm, ierr)
END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE MPI_SYNC(comm)
IMPLICIT NONE
INTEGER :: comm, ierr
CALL MPI_BARRIER(comm, ierr)
END SUBROUTINE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE MPIWaitTurn(comm, myrank, nproc, lPass)
IMPLICIT NONE
INTEGER :: COMM, nProc, Myrank
LOGICAL :: lPass
INTEGER :: i, j, ierr, ito, ifr, istat

IF(.NOT. lPass) THEN
  CALL MPI_BARRIER(Comm, ierr)
  i=0
  DO j = 0, myrank
    CALL MPI_BCAST(i, 1, MPI_INTEGER, j, comm, ierr)
  ENDDO
  !IF(myrank .NE. 0) THEN
  !  ifr = myrank - 1  
    !CALL MPI_RECV(i, 1, MPI_INTEGER, ifr, 0, comm, istat, ierr)
  !ENDIF
ELSE
  DO j = myrank+1, nproc-1
    CALL MPI_BCAST(i, 1, MPI_INTEGER, j, comm, ierr)
  ENDDO
  !IF(myrank .NE. nproc-1) THEN
  !  i = 1
  !  ito = myrank + 1
  !  CALL MPI_SEND(i, 1, MPI_INTEGER, ito, 0, comm, ierr)
  !ENDIF
  CALL MPI_BARRIER(Comm, ierr)
ENDIF
END SUBROUTINE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE MPI_MAX_REAL(dat, comm, lall)
REAL(8) :: dat
INTEGER :: comm
LOGICAL :: lall
INTEGER :: ierr

REAL(8) :: buf

IF(.NOT. lAll) THEN
  CALL MPI_REDUCE(dat, Buf, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, comm, ierr)
ELSE
  CALL MPI_ALLREDUCE(dat, Buf, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm, ierr)
ENDIF
dat = buf
END SUBROUTINE

SUBROUTINE MPI_MAX_INT(dat, comm, lall)
INTEGER :: dat
INTEGER :: comm
LOGICAL :: lall
INTEGER :: ierr

INTEGER :: buf

IF(.NOT. lAll) THEN
  CALL MPI_REDUCE(dat, Buf, 1, MPI_INTEGER, MPI_MAX, 0, comm, ierr)
ELSE
  CALL MPI_ALLREDUCE(dat, Buf, 1, MPI_INTEGER, MPI_MAX, comm, ierr)
ENDIF
dat = buf
END SUBROUTINE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE REDUCEnBCAST_Real0D(dat, Comm1, Comm2, lMember1, lMember2, lall)
IMPLICIT NONE
REAL(8) :: dat
INTEGER :: comm1, comm2
LOGICAL :: lMember1, lMember2, lall
REAL(8) :: temp
!!!Master of comm1 and Master of Comm2 should be identical
IF(lMember1) THEN
  CALL REDUCE(dat, temp, comm1, lall)
ENDIF

IF(lMember2) THEN
  CALL BCAST(temp, comm2)
  dat = temp
ENDIF
END SUBROUTINE

SUBROUTINE REDUCEnBCAST_INT0D(dat, Comm1, Comm2, lMember1, lMember2, lall)
IMPLICIT NONE
INTEGER :: dat
INTEGER :: comm1, comm2
LOGICAL :: lMember1, lMember2, lall
INTEGER :: temp
!!!Master of comm1 and Master of Comm2 should be identical
IF(lMember1) THEN
  CALL REDUCE(dat, temp, comm1, lall)
ENDIF

IF(lMember2) THEN
  CALL BCAST(temp, comm2)
  dat = temp
ENDIF

END SUBROUTINE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE

#endif
